Calorie Data

Data Preparation

Let’s start by loading in our packages:

library(tidyverse)
library(readxl)
library(rstatix)
library(ggpubr)
library(ez)
library(emmeans)
library(car)
library(Hmisc)
library(corrplot)
library(ggsignif)
library(writexl)

And let’s read our excel data in

#Vector for column names and diet names
column_names <- c("rat", "zero", "twenty", "forty", "sixty", "eighty", "hundred", "h2o", "spill", "condition", "day")
diet_names <- c("zero", "twenty", "forty", "sixty", "eighty", "hundred")
condition_names <- c("HI", "LO", "VEH", "NAL")

#Read in data
data <- read_excel("ID Test Data.xlsx", sheet = 6, range = "A3:K50", col_names = column_names)

#turn negatives to zeroes
data_filtered <- data
data_filtered[data_filtered < 0] <- 0

#filter out ID11 (because we missed on his surgery)
data_filtered <- data_filtered %>%
  filter(rat != "ID11")

Now we need to pivot the data into long format:

data_long <- data_filtered %>%
  pivot_longer(cols = diet_names, names_to = "diet", values_to = "kcal")

Also, because spillage and water consumption are in grams, for now I’m just going to drop it:

data_long_foodOnly <- data_long %>%
  subset(select = -c(h2o, spill))

And with that I think we have the data looking how we want it.

Basic visualizations

Now lets try visualizing things before we do any inferential analyses. For a start, lets plot average caloric intake over drug condition

#mean calories for each drug condition (collapsed across days, conditions, and rats)
meanKcals_groupedByDrug <- data_long_foodOnly %>%
  group_by(condition) %>%
  summarise(kcal = mean(kcal))

#SD's
SDKcals_groupedByDrug <- data_long_foodOnly %>%
  group_by(condition) %>%
  summarise(sd = sd(kcal),
            n = n(),
            se = sd/sqrt(n))

#Merge data frames
data_groupedByDrug <- inner_join(meanKcals_groupedByDrug,
                                 SDKcals_groupedByDrug, 
                                by = "condition" )

data_groupedByDrug$condition <- factor(data_groupedByDrug$condition,
                                       levels = c("HI", "LO", "VEH", "NAL"))

#Make plot
my_drug_labels <- setNames(c("High DAMGO", "Low DAMGO", "Naltrexone", "Saline"),
                        c("HI", "LO", "NAL", "VEH"))

data_groupedByDrug %>%
  ggplot(mapping = aes(x = condition,
                       y = kcal,
                       fill = condition,
                       ymin = kcal - se,
                       ymax = kcal + se)) +
  geom_col(show.legend = FALSE) + 
  geom_errorbar(position = position_dodge(0.9),
                width = 0.5,
                size = 0.2) +
  scale_x_discrete(labels = my_drug_labels) +
  labs(x = "Drug Condition",
       y = "Intake (kcal)") +
  theme_classic()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

And to get a sense of between-subject variability, lets do a similar analysis grouping by rat

#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanKcals_groupedByRat <- data_long_foodOnly %>%
  group_by(rat) %>%
  summarise(mean(kcal))

#SD's
SDKcals_groupedByRat <- data_long_foodOnly %>%
  group_by(rat) %>%
  summarise(sd(kcal))

#Merge data frames
data_groupedByRat <- inner_join(meanKcals_groupedByRat,
                                 SDKcals_groupedByRat, 
                                by = "rat" ) %>%
  rename("mean_kcal" = "mean(kcal)",
         "sd_kcal" = "sd(kcal)")

rat_names <- c("ID1", "ID2", "ID3", "ID4", "ID5", "ID6", "ID7", "ID8", "ID9", "ID10", "ID12")
data_groupedByRat$rat <- factor(data_groupedByRat$rat,
                                       levels = rat_names)

#Make plot
data_groupedByRat %>%
  ggplot(mapping = aes(x = rat,
                       y = mean_kcal,
                       fill = rat)) +
  geom_col() + 
  geom_errorbar(aes(ymin = mean_kcal - sd_kcal,
                    ymax = mean_kcal + sd_kcal))

And just as a rough replication of the last study, lets group by diet too:

#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanKcals_groupedByDiet <- data_long_foodOnly %>%
  group_by(diet) %>%
  summarise(kcal = mean(kcal))

#SD's
SEKcals_groupedByDiet <- data_long_foodOnly %>%
  group_by(diet) %>%
  summarise(sd = sd(kcal),
            n = n(),
            se = sd/sqrt(n))

#Merge data frames
data_groupedByDiet <- inner_join(meanKcals_groupedByDiet,
                                 SEKcals_groupedByDiet, 
                                by = "diet" )

data_groupedByDiet$diet <- factor(data_groupedByDiet$diet,
                                       levels = diet_names)

#Make plot
my_diet_labels <- setNames(c("0%", "20%", "40%", "60%", "80%", "100%"),
                           diet_names)


data_groupedByDiet %>%
  ggplot(mapping = aes(x = diet,
                       y = kcal,
                       fill = diet,
                       ymin = kcal - se,
                       ymax = kcal + se)) +
  geom_col(show.legend = FALSE) + 
  geom_errorbar(position = position_dodge(0.9),
                width = 0.5,
                size = 0.2) +
  scale_x_discrete(labels = my_diet_labels) +
  labs(x = "Diet (% sugar by weight)",
       y = "Intake (kcal)") +
  theme_classic()

So here it looks like on average the eighty percent diet was a hit. However, it seems like across all of these dimensions, as overall caloric intake increases, the variation in intake also increases.

Inferential Stats

Now let’s do the hypothesis test:

res.aov = anova_test(
  data = data_long_foodOnly, dv = kcal, wid = rat,
  within = c(condition, diet)
)
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##           Effect   DFn    DFd      F        p p<.05   ges
## 1      condition  3.00  30.00 29.188 4.95e-09     * 0.070
## 2           diet  2.17  21.71  7.929 2.00e-03     * 0.313
## 3 condition:diet 15.00 150.00  4.536 4.66e-07     * 0.145
#trying it a different way
exp2 <- aov(kcal ~ diet*condition + Error(rat/(diet*condition)), data = data_long_foodOnly)
summary(exp2)
## 
## Error: rat
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 10  980.2   98.02               
## 
## Error: rat:diet
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## diet       5  17228    3446   7.929 1.46e-05 ***
## Residuals 50  21727     435                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: rat:condition
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## condition  3 2855.2   951.7   29.19 4.95e-09 ***
## Residuals 30  978.2    32.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: rat:diet:condition
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## diet:condition  15   6394   426.2   4.536 4.66e-07 ***
## Residuals      150  14094    94.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#And a third way because the first two don't match
kcals.aov <- ezANOVA(data = data_long_foodOnly,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(diet, condition),
                    return_aov = TRUE,
                    type = 1)
## Warning: Converting "rat" to factor for ANOVA.
## Warning: Converting "diet" to factor for ANOVA.
## Warning: Converting "condition" to factor for ANOVA.
kcals.aov$ANOVA
##           Effect DFn DFd         F            p p<.05        ges
## 1           diet   5  50  7.929407 1.464019e-05     * 0.31887886
## 2      condition   3  30 29.187765 4.948465e-09     * 0.07200046
## 3 diet:condition  15 150  4.536271 4.661740e-07     * 0.14802151

So these two methods of calculation slightly disagree about the p-value for the main effect of drug, but they both clearly show two main effects and an interaction. The p-values are actually suspiciously small here; we might have to circle back to this.

In any case, I want to get a more detailed visualization before I decide what post-hoc analyses I want to run:

ggboxplot(data_long_foodOnly,
          x = "diet", y = "kcal",
          color = "condition"
          )

#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanKcals_groupedByDiet <- data_long_foodOnly %>%
  group_by(diet, condition) %>%
  summarise(kcal = mean(kcal))
## `summarise()` has grouped output by 'diet'. You can override using the
## `.groups` argument.
#SD's
SEKcals_groupedByDiet <- data_long_foodOnly %>%
  group_by(diet, condition) %>%
  summarise(n = n(),
            sd = sd(kcal),
            se = sd/sqrt(n))
## `summarise()` has grouped output by 'diet'. You can override using the
## `.groups` argument.
#Merge data frames
data_groupedByDiet <- inner_join(meanKcals_groupedByDiet,
                                 SEKcals_groupedByDiet, 
                                 by = c("diet", "condition"))
  

data_groupedByDiet$diet <- factor(data_groupedByDiet$diet,
                                       levels = diet_names)

#Make plot
my_x_labels <- setNames(c("High DAMGO", "Low DAMGO", "Naltrexone", "Saline"),
                        c("HI", "LO", "NAL", "VEH"))
my_fill_labels <- setNames(c("0%", "20%", "40%", "60%", "80%", "100%"),
                           diet_names)



data_groupedByDiet %>%
  ggplot(mapping = aes(x = condition,
                       y = kcal,
                       fill = diet,
                       ymin = kcal - se,
                       ymax = kcal + se)) +
  geom_col(position = "dodge") + 
  geom_errorbar(position = position_dodge(0.9),
                width = 0.5,
                size = 0.2) + 
  scale_x_discrete(labels = my_x_labels) +
  scale_fill_discrete(labels = my_fill_labels) +
  labs(x = "Drug Condition",
       y = "Intake (kcal)",
       fill = "Diet (% sugar)") +
  theme_classic() +
  theme_classic(base_size = 16) +
  theme(axis.text.x = element_text(size = 12))

So at least descriptively, this looks like good news! Consumption was highest of the sixty and eighty percent diet, and within each group, drugs look like they behaved the way we were expecting. Additionally, the significant interaction suggests that the drug effect was diet-dependent; the high dose of DAMGO spiked calorie intake a lot more for the eighty percent diet than it did for any of the others.

However, what I’m concerned about is that calculating things based on calories might actually be giving us an underestimate of what’s going on. That’s because the diets favored by the rats are actually the ones with lower caloric density – the ones with higher sugar percentages by weight. I want to check this by repeating the same analysis pipeline from the data in grams.

Weight Data

Data preparation

#Vector for column names and diet names
column_names <- c("rat", "zero", "twenty", "forty", "sixty", "eighty", "hundred", "h2o", "spill", "condition", "day")
diet_names <- c("zero", "twenty", "forty", "sixty", "eighty", "hundred")

#Read in data
data_grams <- read_excel("ID Test Data.xlsx", sheet = 5, range = "A3:K50", col_names = column_names)

#turn negatives to zeroes
data_grams_filtered <- data_grams
data_grams_filtered[data_grams_filtered < 0] <- 0

#filter out ID11
data_grams_filtered <-data_grams_filtered %>%
  filter(rat != "ID11")

#Pivot long
data_grams_long <- data_grams_filtered %>%
  pivot_longer(cols = diet_names, names_to = "diet", values_to = "grams")

#Drop spill and water
data_grams_long_foodOnly <- data_grams_long %>%
  subset(select = -c(h2o, spill))

Basic visualizations

#mean calories for each drug condition (collapsed across days, conditions, and rats)
meanGrams_groupedByDrug <- data_grams_long_foodOnly %>%
  group_by(condition) %>%
  summarise(grams = mean(grams))

#SD's
SDGrams_groupedByDrug <- data_grams_long_foodOnly %>%
  group_by(condition) %>%
  summarise(sd = sd(grams),
            n = n(),
            se = sd/sqrt(n))

#Merge data frames
data_grams_groupedByDrug <- inner_join(meanGrams_groupedByDrug,
                                 SDGrams_groupedByDrug, 
                                by = "condition" )

data_grams_groupedByDrug$condition <- factor(data_grams_groupedByDrug$condition,
                                       levels = c("HI", "LO", "VEH", "NAL"))

#Make plot
data_grams_groupedByDrug %>%
  ggplot(mapping = aes(x = condition,
                       y = grams,
                       fill = condition,
                       ymin = grams - se,
                       ymax = grams + se)) +
  geom_col(show.legend = FALSE) + 
  geom_errorbar(position = position_dodge(0.9),
                width = 0.5,
                size = 0.2) +
  scale_x_discrete(labels = my_drug_labels) +
  labs(x = "Drug Condition",
       y = "Intake (grams)") +
  theme_classic()

#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanGrams_groupedByRat <- data_grams_long_foodOnly %>%
  group_by(rat) %>%
  summarise(mean(grams))

#SD's
SDGrams_groupedByRat <- data_grams_long_foodOnly %>%
  group_by(rat) %>%
  summarise(sd(grams))

#Merge data frames
data_grams_groupedByRat <- inner_join(meanGrams_groupedByRat,
                                 SDGrams_groupedByRat, 
                                by = "rat" ) %>%
  rename("mean_grams" = "mean(grams)",
         "sd_grams" = "sd(grams)")

rat_names <- c("ID1", "ID2", "ID3", "ID4", "ID5", "ID6", "ID7", "ID8", "ID9", "ID10", "ID11", "ID12")
data_grams_groupedByRat$rat <- factor(data_grams_groupedByRat$rat,
                                       levels = rat_names)

#Make plot
data_grams_groupedByRat %>%
  ggplot(mapping = aes(x = rat,
                       y = mean_grams,
                       fill = rat)) +
  geom_col() + 
  geom_errorbar(aes(ymin = mean_grams - sd_grams,
                    ymax = mean_grams + sd_grams))

#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanGrams_groupedByDiet <- data_grams_long_foodOnly %>%
  group_by(diet) %>%
  summarise(grams = mean(grams))

#SD's
SDGrams_groupedByDiet <- data_grams_long_foodOnly %>%
  group_by(diet) %>%
  summarise(sd = sd(grams),
            n = n(),
            se = sd/sqrt(n))

#Merge data frames
data_grams_groupedByDiet <- inner_join(meanGrams_groupedByDiet,
                                 SDGrams_groupedByDiet, 
                                by = "diet" ) 

data_grams_groupedByDiet$diet <- factor(data_grams_groupedByDiet$diet,
                                       levels = diet_names)

#Make plot
data_grams_groupedByDiet %>%
  ggplot(mapping = aes(x = diet,
                       y = grams,
                       fill = diet,
                       ymin = grams - se,
                       ymax = grams + se)) +
  geom_col(show.legend = FALSE) + 
  geom_errorbar(position = position_dodge(0.9),
                width = 0.5,
                size = 0.2) +
  scale_x_discrete(labels = my_diet_labels) +
  labs(x = "Diet (% sugar by weight)",
       y = "Intake (grams)") +
  theme_classic()

Inferential Stats

res.aovGrams = anova_test(
  data = data_grams_long_foodOnly, dv = grams, wid = rat,
  within = c(condition, diet)
)
get_anova_table(res.aovGrams)
## ANOVA Table (type III tests)
## 
##           Effect   DFn    DFd      F        p p<.05   ges
## 1      condition  3.00  30.00 30.621 2.90e-09     * 0.080
## 2           diet  1.83  18.32 10.690 1.00e-03     * 0.367
## 3 condition:diet 15.00 150.00  5.075 4.87e-08     * 0.171
#trying it a different way
exp2Grams <- aov(grams ~ diet*condition + Error(rat/(diet*condition)), data = data_grams_long_foodOnly)
summary(exp2Grams) 
## 
## Error: rat
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 10  28.04   2.804               
## 
## Error: rat:diet
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## diet       5  698.6  139.72   10.69 5.06e-07 ***
## Residuals 50  653.5   13.07                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: rat:condition
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## condition  3 104.39   34.80   30.62 2.9e-09 ***
## Residuals 30  34.09    1.14                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: rat:diet:condition
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## diet:condition  15  247.3  16.490   5.075 4.87e-08 ***
## Residuals      150  487.4   3.249                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggboxplot(data_grams_long_foodOnly,
          x = "diet", y = "grams",
          color = "condition"
          )

#And a third way because the first two don't match
grams.aov <- ezANOVA(data = data_grams_long_foodOnly,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(diet, condition),
                    return_aov = TRUE,
                    type = 1)
## Warning: Converting "rat" to factor for ANOVA.
## Warning: Converting "diet" to factor for ANOVA.
## Warning: Converting "condition" to factor for ANOVA.
grams.aov$ANOVA
##           Effect DFn DFd         F            p p<.05       ges
## 1           diet   5  50 10.690176 5.057562e-07     * 0.3728687
## 2      condition   3  30 30.621458 2.902644e-09     * 0.0815962
## 3 diet:condition  15 150  5.074863 4.872602e-08     * 0.1739020
#Figure 2

#mean calories for each rat (collapsed across days, conditions, and drug condition)
meanGrams <- data_grams_long_foodOnly %>%
  group_by(diet, condition) %>%
  summarise(grams = mean(grams))
## `summarise()` has grouped output by 'diet'. You can override using the
## `.groups` argument.
#SD's
SEGrams <- data_grams_long_foodOnly %>%
  group_by(diet, condition) %>%
  summarise(n = n(),
            sd = sd(grams),
            se = sd/sqrt(n))
## `summarise()` has grouped output by 'diet'. You can override using the
## `.groups` argument.
#Merge data frames
data_joined <- inner_join(meanGrams,
                                 SEGrams, 
                                 by = c("diet", "condition"))
  

data_joined$diet <- factor(data_joined$diet,
                                       levels = diet_names)

#Make plot
my_x_labels <- setNames(c("High DAMGO", "Low DAMGO", "Naltrexone", "Saline"),
                        c("HI", "LO", "NAL", "VEH"))
my_fill_labels <- setNames(c("0%", "20%", "40%", "60%", "80%", "100%"),
                           diet_names)

data_joined %>%
  ggplot(mapping = aes(x = condition,
                       y = grams,
                       fill = diet,
                       ymin = grams - se,
                       ymax = grams+ se)) +
  geom_col(position = "dodge") + 
  geom_errorbar(position = position_dodge(0.9),
                width = 0.5,
                size = 0.2) + 
  scale_x_discrete(labels = my_x_labels) +
  scale_fill_discrete(labels = my_fill_labels) +
  labs(x = "Drug Condition",
       y = "Intake (grams)",
       fill = "Diet (% sugar)") +
  theme_classic(base_size = 16) +
  theme(axis.text.x = element_text(size = 12))

So the story in terms of grams looks just about the same as it does in terms of calories. However, there’s enough variability in the data that two possibilities are jumping into my head:

1.) Each individual rat may have very different preferences. Some may be sugar-preferrers while others are fat-preferrers 2.) Outliers may be driving some of the results we’re seeing. For instance, we had one rat, eat 23 grams of the 80% diet on one day, and another ate 17. I feel like we need to account for that somehow.

Before I try to tackle this, however, let’s get into some post-hoc analyses.

Post Hoc Analyses

The challenge throughout this section is going to be to capture what I want to see without running a million tests. To a certain extent I’m not sure that’ll be possible, but I do think I should make an effort.

Looking back at my prospectus, the questions I want to answer are:

1.) What ratio do the rats prefer (should be a replication from experiment 1)

2.) Do the rats prefer blend diets to the single nutrient ones (again, a replication)

3.) Does opioid blockade shift diet preference to lower fat or abolish preferences altogether?

4.) Does opioid stimulation shift preference to higher fat?

5.) Do opioids overall increase consumption?

6.) Does naltrexone overall decrease consumption?

Question 1

To figure out which diets the rats preferred overall, let’s first visualize the overall pattern. We actually did a similar plot earlier, but lets try a boxplot this time:

#For calories
data_long_foodOnly %>%
ggplot(mapping = aes(x = factor(diet, level = diet_names), 
                     y = kcal, fill = diet)) + 
         geom_boxplot()

#For grams
data_grams_long_foodOnly %>%
ggplot(mapping = aes(x = factor(diet, level = diet_names), 
                     y = grams, fill = diet)) + 
         geom_boxplot()

And now that we have a visual, we can run pairwise paired t-tests between diet conditions. To be super conservative, let’s start with Bonferroni corrections

#By calories
pairwiseKcal <- data_long_foodOnly %>%
  pairwise_t_test(kcal~diet, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
pairwiseKcal
## # A tibble: 15 × 10
##    .y.   group1  group2    n1    n2 statistic    df       p   p.adj p.adj.signif
##  * <chr> <chr>   <chr>  <int> <int>     <dbl> <dbl>   <dbl>   <dbl> <chr>       
##  1 kcal  eighty  forty     44    44    5.66      43 1.13e-6 1.7 e-5 ****        
##  2 kcal  eighty  hundr…    44    44    6.01      43 3.52e-7 5.28e-6 ****        
##  3 kcal  eighty  sixty     44    44    1.82      43 7.5 e-2 1   e+0 ns          
##  4 kcal  eighty  twenty    44    44    5.00      43 1.02e-5 1.53e-4 ***         
##  5 kcal  eighty  zero      44    44    5.20      43 5.27e-6 7.91e-5 ****        
##  6 kcal  forty   hundr…    44    44    1.10      43 2.76e-1 1   e+0 ns          
##  7 kcal  forty   sixty     44    44   -3.96      43 2.79e-4 4   e-3 **          
##  8 kcal  forty   twenty    44    44   -1.05      43 2.99e-1 1   e+0 ns          
##  9 kcal  forty   zero      44    44   -0.986     43 3.3 e-1 1   e+0 ns          
## 10 kcal  hundred sixty     44    44   -4.20      43 1.31e-4 2   e-3 **          
## 11 kcal  hundred twenty    44    44   -1.45      43 1.55e-1 1   e+0 ns          
## 12 kcal  hundred zero      44    44   -1.53      43 1.34e-1 1   e+0 ns          
## 13 kcal  sixty   twenty    44    44    3.27      43 2   e-3 3.2 e-2 *           
## 14 kcal  sixty   zero      44    44    3.15      43 3   e-3 4.4 e-2 *           
## 15 kcal  twenty  zero      44    44    0.0510    43 9.6 e-1 1   e+0 ns
#By grams
pairwiseGrams <- data_grams_long_foodOnly %>%
  pairwise_t_test(grams~diet, paired = TRUE,
                  p.adjust.method = "bonferroni")
pairwiseGrams
## # A tibble: 15 × 10
##    .y.   group1  group2    n1    n2 statistic    df       p   p.adj p.adj.signif
##  * <chr> <chr>   <chr>  <int> <int>     <dbl> <dbl>   <dbl>   <dbl> <chr>       
##  1 grams eighty  forty     44    44    5.90      43 5.16e-7 7.74e-6 ****        
##  2 grams eighty  hundr…    44    44    5.88      43 5.46e-7 8.19e-6 ****        
##  3 grams eighty  sixty     44    44    2.42      43 2   e-2 2.96e-1 ns          
##  4 grams eighty  twenty    44    44    5.57      43 1.52e-6 2.28e-5 ****        
##  5 grams eighty  zero      44    44    5.80      43 7.25e-7 1.09e-5 ****        
##  6 grams forty   hundr…    44    44   -0.0743    43 9.41e-1 1   e+0 ns          
##  7 grams forty   sixty     44    44   -4.08      43 1.9 e-4 3   e-3 **          
##  8 grams forty   twenty    44    44   -0.754     43 4.55e-1 1   e+0 ns          
##  9 grams forty   zero      44    44   -0.357     43 7.23e-1 1   e+0 ns          
## 10 grams hundred sixty     44    44   -3.88      43 3.57e-4 5   e-3 **          
## 11 grams hundred twenty    44    44   -0.401     43 6.9 e-1 1   e+0 ns          
## 12 grams hundred zero      44    44   -0.156     43 8.77e-1 1   e+0 ns          
## 13 grams sixty   twenty    44    44    3.67      43 6.62e-4 1   e-2 **          
## 14 grams sixty   zero      44    44    3.68      43 6.4 e-4 1   e-2 **          
## 15 grams twenty  zero      44    44    0.283     43 7.79e-1 1   e+0 ns

So even being as conservative as possible (which we may decide is what we want) the 80% diet is significantly more consumed than all other diets except the 60%, and the 60% diet is significantly preferred over the 0%, 40%, and 100% diets.There apparently isn’t a significant difference between the 20% and 60% diets, however, which actually makes me want to take a look at the effect sizes we’re dealing with:

#By calories
data_long_foodOnly %>%
  cohens_d(kcal~diet, paired = TRUE)
## # A tibble: 15 × 7
##    .y.   group1  group2   effsize    n1    n2 magnitude 
##  * <chr> <chr>   <chr>      <dbl> <int> <int> <ord>     
##  1 kcal  eighty  forty    0.854      44    44 large     
##  2 kcal  eighty  hundred  0.906      44    44 large     
##  3 kcal  eighty  sixty    0.275      44    44 small     
##  4 kcal  eighty  twenty   0.753      44    44 moderate  
##  5 kcal  eighty  zero     0.784      44    44 moderate  
##  6 kcal  forty   hundred  0.166      44    44 negligible
##  7 kcal  forty   sixty   -0.597      44    44 moderate  
##  8 kcal  forty   twenty  -0.158      44    44 negligible
##  9 kcal  forty   zero    -0.149      44    44 negligible
## 10 kcal  hundred sixty   -0.633      44    44 moderate  
## 11 kcal  hundred twenty  -0.218      44    44 small     
## 12 kcal  hundred zero    -0.230      44    44 small     
## 13 kcal  sixty   twenty   0.493      44    44 small     
## 14 kcal  sixty   zero     0.475      44    44 small     
## 15 kcal  twenty  zero     0.00768    44    44 negligible
#By grams
data_grams_long_foodOnly %>%
  cohens_d(grams~diet, paired = TRUE)
## # A tibble: 15 × 7
##    .y.   group1  group2  effsize    n1    n2 magnitude 
##  * <chr> <chr>   <chr>     <dbl> <int> <int> <ord>     
##  1 grams eighty  forty    0.889     44    44 large     
##  2 grams eighty  hundred  0.886     44    44 large     
##  3 grams eighty  sixty    0.365     44    44 small     
##  4 grams eighty  twenty   0.840     44    44 large     
##  5 grams eighty  zero     0.874     44    44 large     
##  6 grams forty   hundred -0.0112    44    44 negligible
##  7 grams forty   sixty   -0.615     44    44 moderate  
##  8 grams forty   twenty  -0.114     44    44 negligible
##  9 grams forty   zero    -0.0538    44    44 negligible
## 10 grams hundred sixty   -0.584     44    44 moderate  
## 11 grams hundred twenty  -0.0605    44    44 negligible
## 12 grams hundred zero    -0.0236    44    44 negligible
## 13 grams sixty   twenty   0.553     44    44 moderate  
## 14 grams sixty   zero     0.555     44    44 moderate  
## 15 grams twenty  zero     0.0426    44    44 negligible

So they’re medium to large, which is encouraging to see. Based on these analyses, I feel pretty comfortable concluding that the 80% diet was the favorite on average.

Question 2

Based on the plots we’ve generated, it’s pretty clear that the rats preferred blend diets on average over single-nutrient ones. However, just to be above reproach I want to do a formal comparison. Time for fun with contrasts! Assigning a coefficient of 1/4 for all the blend diets and a coefficient of -1/2 for the single diet conditions should give us what we want.

First I have to do some R trickery to force the software to handle my variables the way I want it to:

#kcals
data_long_foodOnly_factor <- data_long_foodOnly

names <- c(1:4)

data_long_foodOnly_factor[,names] <- lapply(data_long_foodOnly_factor[,names], factor)

data_long_foodOnly_factor$rat <- factor(data_long_foodOnly_factor$rat, levels = rat_names)

data_long_foodOnly_factor$diet <- factor(data_long_foodOnly_factor$diet, levels = diet_names)

data_long_foodOnly_factor$condition <- factor(data_long_foodOnly_factor$condition, levels = condition_names)

str(data_long_foodOnly_factor)
## tibble [264 × 5] (S3: tbl_df/tbl/data.frame)
##  $ rat      : Factor w/ 12 levels "ID1","ID2","ID3",..: 1 1 1 1 1 1 2 2 2 2 ...
##  $ condition: Factor w/ 4 levels "HI","LO","VEH",..: 4 4 4 4 4 4 2 2 2 2 ...
##  $ day      : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ diet     : Factor w/ 6 levels "zero","twenty",..: 1 2 3 4 5 6 1 2 3 4 ...
##  $ kcal     : num [1:264] 0.45 1.76 3.92 6.3 10.35 ...
#grams
data_grams_long_foodOnly_factor <- data_grams_long_foodOnly

names <- c(1:4)

data_grams_long_foodOnly_factor[,names] <- lapply(data_grams_long_foodOnly_factor[,names], factor)

data_grams_long_foodOnly_factor$rat <- factor(data_grams_long_foodOnly_factor$rat, levels = rat_names)

data_grams_long_foodOnly_factor$diet <- factor(data_grams_long_foodOnly_factor$diet, levels = diet_names)

data_grams_long_foodOnly_factor$condition <- factor(data_grams_long_foodOnly_factor$condition, levels = condition_names)

str(data_grams_long_foodOnly_factor)
## tibble [264 × 5] (S3: tbl_df/tbl/data.frame)
##  $ rat      : Factor w/ 12 levels "ID1","ID2","ID3",..: 1 1 1 1 1 1 2 2 2 2 ...
##  $ condition: Factor w/ 4 levels "HI","LO","VEH",..: 4 4 4 4 4 4 2 2 2 2 ...
##  $ day      : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ diet     : Factor w/ 6 levels "zero","twenty",..: 1 2 3 4 5 6 1 2 3 4 ...
##  $ grams    : num [1:264] 0.05 0.22 0.56 1.05 2.07 ...

And now I’m going to give it a shot:

#Set contrast coefficient matrix
#By kcals
options(contrasts = c("contr.sum", "contr.poly"))
kcal.aov <- ezANOVA(data = data_long_foodOnly_factor,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(diet, condition),
                    return_aov = TRUE,
                    type = 1)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
kcal.em <- emmeans(kcal.aov$aov, ~diet * condition)

contrast_names <- list(c(2, -1, -1, -1, -1, 2,
                         2, -1, -1, -1, -1, 2,
                         2, -1, -1, -1, -1, 2,
                         2, -1, -1, -1, -1, 2))

contrast(kcal.em,
         method = contrast_names)
##  contrast                                                      estimate   SE df
##  c(2, -1, -1, -1, -1, 2, 2, -1, -1, -1, -1, 2, 2, -1, -1, -1,      -134 43.5 50
##  t.ratio p.value
##   -3.083  0.0033
#By grams
options(contrasts = c("contr.sum", "contr.poly"))
grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(diet, condition),
                    return_aov = TRUE,
                    type = 1)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
grams.em <- emmeans(grams.aov$aov, ~diet * condition)

contrast_names <- list(c(2, -1, -1, -1, -1, 2,
                         2, -1, -1, -1, -1, 2,
                         2, -1, -1, -1, -1, 2,
                         2, -1, -1, -1, -1, 2))

contrast(grams.em,
         method = contrast_names)
##  contrast                                                      estimate   SE df
##  c(2, -1, -1, -1, -1, 2, 2, -1, -1, -1, -1, 2, 2, -1, -1, -1,     -25.3 7.55 50
##  t.ratio p.value
##   -3.348  0.0016

So if I did this right, then we did find the effect we were looking for. However, I’m a little unsure of this, so I may double-check in SPSS just for peace of mind.

Question 3

To go about answering this question, we can start by filtering the data to only include the naltrexone and saline groups:

#kcals
data_long_foodOnly_factor_q3 <- data_long_foodOnly_factor %>%
  filter(condition == "VEH" | condition == "NAL")

#grams
data_grams_long_foodOnly_factor_q3 <- data_grams_long_foodOnly_factor %>%
   filter(condition == "VEH" | condition == "NAL")

And now we can visualize:

#kcals
ggboxplot(data_long_foodOnly_factor_q3,
          x = "condition", y = "kcal",
          color = "diet"
          )

#grams
ggboxplot(data_grams_long_foodOnly_factor_q3,
          x = "condition", y = "grams",
          color = "diet"
          )

So just visually, it looks like the trend of consumption rising as sugar content increases holds for both the vehicle and naltrexone groups.

Since we have a significant omnibus effect, I feel comfortable proceeding with looking at significant one-way effects within each drug category. We can start by filtering to create datasets for each specific drug group:

## kcals
# VEH dataset
data_long_foodOnly_factor_veh <- data_long_foodOnly_factor %>%
  filter(condition == "VEH")

#NAL dataset
data_long_foodOnly_factor_nal <- data_long_foodOnly_factor %>%
  filter(condition == "NAL")

##grams
# VEH dataset
data_grams_long_foodOnly_factor_veh <- data_grams_long_foodOnly_factor %>%
  filter(condition == "VEH")

#NAL dataset
data_grams_long_foodOnly_factor_nal <- data_grams_long_foodOnly_factor %>%
  filter(condition == "NAL")

Now we can run one-way ANOVA’s on these data:

## kcals
#VEH
veh.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_veh,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(diet),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
veh.kcal.aov$ANOVA
##   Effect DFn DFd        F            p p<.05      ges
## 2   diet   5  50 8.839876 4.603764e-06     * 0.434701
#NAL
nal.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_nal,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(diet),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
nal.kcal.aov$ANOVA
##   Effect DFn DFd        F          p p<.05       ges
## 2   diet   5  50 2.328814 0.05607767       0.1798096
## grams
veh.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_veh,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(diet),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
veh.grams.aov$ANOVA
##   Effect DFn DFd        F            p p<.05       ges
## 2   diet   5  50 12.29994 8.511631e-08     * 0.5238583
#NAL
nal.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_nal,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(diet),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
nal.grams.aov$ANOVA
##   Effect DFn DFd        F           p p<.05       ges
## 2   diet   5  50 3.677905 0.006531315     * 0.2582462

So we have an omnibus effect of diet on consumption for the vehicle condition, but not the naltrexone group. Within the vehicle condition, we can proceed with pairwise comparisons:

###VEH
#Significance
data_long_foodOnly_factor_veh %>%
  pairwise_t_test(kcal~diet, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
## # A tibble: 15 × 10
##    .y.   group1 group2     n1    n2 statistic    df        p p.adj p.adj.signif
##  * <chr> <chr>  <chr>   <int> <int>     <dbl> <dbl>    <dbl> <dbl> <chr>       
##  1 kcal  zero   twenty     11    11    -0.695    10 0.503    1     ns          
##  2 kcal  zero   forty      11    11    -0.704    10 0.498    1     ns          
##  3 kcal  zero   sixty      11    11    -2.10     10 0.062    0.928 ns          
##  4 kcal  zero   eighty     11    11    -4.29     10 0.002    0.024 *           
##  5 kcal  zero   hundred    11    11     2.06     10 0.066    0.996 ns          
##  6 kcal  twenty forty      11    11     0.470    10 0.648    1     ns          
##  7 kcal  twenty sixty      11    11    -1.21     10 0.253    1     ns          
##  8 kcal  twenty eighty     11    11    -3.97     10 0.003    0.04  *           
##  9 kcal  twenty hundred    11    11     1.34     10 0.21     1     ns          
## 10 kcal  forty  sixty      11    11    -1.90     10 0.087    1     ns          
## 11 kcal  forty  eighty     11    11    -4.07     10 0.002    0.034 *           
## 12 kcal  forty  hundred    11    11     2.64     10 0.025    0.37  ns          
## 13 kcal  sixty  eighty     11    11    -1.69     10 0.122    1     ns          
## 14 kcal  sixty  hundred    11    11     2.78     10 0.019    0.29  ns          
## 15 kcal  eighty hundred    11    11     5.26     10 0.000366 0.005 **
#Effect size
data_long_foodOnly_factor_veh %>%
  cohens_d(kcal~diet, paired = TRUE)
## # A tibble: 15 × 7
##    .y.   group1 group2  effsize    n1    n2 magnitude 
##  * <chr> <chr>  <chr>     <dbl> <int> <int> <ord>     
##  1 kcal  zero   twenty   -0.209    11    11 small     
##  2 kcal  zero   forty    -0.212    11    11 small     
##  3 kcal  zero   sixty    -0.634    11    11 moderate  
##  4 kcal  zero   eighty   -1.29     11    11 large     
##  5 kcal  zero   hundred   0.621    11    11 moderate  
##  6 kcal  twenty forty     0.142    11    11 negligible
##  7 kcal  twenty sixty    -0.366    11    11 small     
##  8 kcal  twenty eighty   -1.20     11    11 large     
##  9 kcal  twenty hundred   0.404    11    11 small     
## 10 kcal  forty  sixty    -0.572    11    11 moderate  
## 11 kcal  forty  eighty   -1.23     11    11 large     
## 12 kcal  forty  hundred   0.796    11    11 moderate  
## 13 kcal  sixty  eighty   -0.509    11    11 moderate  
## 14 kcal  sixty  hundred   0.839    11    11 large     
## 15 kcal  eighty hundred   1.59     11    11 large
###NAL
##kcals
#Significance
data_long_foodOnly_factor_nal %>%
  pairwise_t_test(kcal~diet, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
## # A tibble: 15 × 10
##    .y.   group1 group2     n1    n2 statistic    df     p p.adj p.adj.signif
##  * <chr> <chr>  <chr>   <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
##  1 kcal  zero   twenty     11    11     0.249    10 0.808 1     ns          
##  2 kcal  zero   forty      11    11     1.13     10 0.287 1     ns          
##  3 kcal  zero   sixty      11    11    -0.782    10 0.453 1     ns          
##  4 kcal  zero   eighty     11    11    -1.59     10 0.143 1     ns          
##  5 kcal  zero   hundred    11    11     1.46     10 0.174 1     ns          
##  6 kcal  twenty forty      11    11     1.01     10 0.337 1     ns          
##  7 kcal  twenty sixty      11    11    -0.977    10 0.352 1     ns          
##  8 kcal  twenty eighty     11    11    -2.04     10 0.069 1     ns          
##  9 kcal  twenty hundred    11    11     1.21     10 0.253 1     ns          
## 10 kcal  forty  sixty      11    11    -1.60     10 0.142 1     ns          
## 11 kcal  forty  eighty     11    11    -3.89     10 0.003 0.045 *           
## 12 kcal  forty  hundred    11    11     1.58     10 0.146 1     ns          
## 13 kcal  sixty  eighty     11    11    -0.340    10 0.741 1     ns          
## 14 kcal  sixty  hundred    11    11     1.81     10 0.101 1     ns          
## 15 kcal  eighty hundred    11    11     4.47     10 0.001 0.018 *
#Effect size
data_long_foodOnly_factor_nal %>%
  cohens_d(kcal~diet, paired = TRUE)
## # A tibble: 15 × 7
##    .y.   group1 group2  effsize    n1    n2 magnitude 
##  * <chr> <chr>  <chr>     <dbl> <int> <int> <ord>     
##  1 kcal  zero   twenty   0.0751    11    11 negligible
##  2 kcal  zero   forty    0.339     11    11 small     
##  3 kcal  zero   sixty   -0.236     11    11 small     
##  4 kcal  zero   eighty  -0.479     11    11 small     
##  5 kcal  zero   hundred  0.441     11    11 small     
##  6 kcal  twenty forty    0.304     11    11 small     
##  7 kcal  twenty sixty   -0.295     11    11 small     
##  8 kcal  twenty eighty  -0.614     11    11 moderate  
##  9 kcal  twenty hundred  0.366     11    11 small     
## 10 kcal  forty  sixty   -0.481     11    11 small     
## 11 kcal  forty  eighty  -1.17      11    11 large     
## 12 kcal  forty  hundred  0.475     11    11 small     
## 13 kcal  sixty  eighty  -0.103     11    11 negligible
## 14 kcal  sixty  hundred  0.545     11    11 moderate  
## 15 kcal  eighty hundred  1.35      11    11 large
##grams
#Significance
data_grams_long_foodOnly_factor_nal %>%
  pairwise_t_test(grams~diet, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
## # A tibble: 15 × 10
##    .y.   group1 group2     n1    n2 statistic    df     p p.adj p.adj.signif
##  * <chr> <chr>  <chr>   <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
##  1 grams zero   twenty     11    11     0.125    10 0.903 1     ns          
##  2 grams zero   forty      11    11     0.994    10 0.344 1     ns          
##  3 grams zero   sixty      11    11    -1.12     10 0.29  1     ns          
##  4 grams zero   eighty     11    11    -2.64     10 0.025 0.374 ns          
##  5 grams zero   hundred    11    11     1.33     10 0.213 1     ns          
##  6 grams twenty forty      11    11     0.947    10 0.366 1     ns          
##  7 grams twenty sixty      11    11    -1.19     10 0.26  1     ns          
##  8 grams twenty eighty     11    11    -2.89     10 0.016 0.243 ns          
##  9 grams twenty hundred    11    11     1.07     10 0.31  1     ns          
## 10 grams forty  sixty      11    11    -1.64     10 0.132 1     ns          
## 11 grams forty  eighty     11    11    -4.11     10 0.002 0.032 *           
## 12 grams forty  hundred    11    11     1.03     10 0.329 1     ns          
## 13 grams sixty  eighty     11    11    -0.666    10 0.521 1     ns          
## 14 grams sixty  hundred    11    11     1.76     10 0.109 1     ns          
## 15 grams eighty hundred    11    11     4.42     10 0.001 0.019 *

So overall, it appears that the sugarier diets are favored in the vehicle group, but these preferences are abolished (or at least are drastically reduced) in the naltrexone condition.

Question 4

I think I can repeat this pipeline to compare what happens between the DAMGO groups and the vehicle group.

Always a good idea to start with visualization:

# kcals
data_long_foodOnly_factor_q3 <- data_long_foodOnly_factor %>%
  filter(condition == "HI" | condition == "LO" | condition == "VEH")

ggboxplot(data_long_foodOnly_factor_q3,
          x = "condition", y = "kcal",
          color = "diet"
          )

# grams
data_grams_long_foodOnly_factor_q3 <- data_grams_long_foodOnly_factor %>%
  filter(condition == "HI" | condition == "LO" | condition == "VEH")

ggboxplot(data_grams_long_foodOnly_factor_q3,
          x = "condition", y = "grams",
          color = "diet"
          )

Again, the significant omnibus effect allows us to run simple one-way ANOVAs for the HI and LO dose groups:

###Filter out HI and LO groups
#kcals
# HI dataset
data_long_foodOnly_factor_hi <- data_long_foodOnly_factor %>%
  filter(condition == "HI")

#LO dataset
data_long_foodOnly_factor_lo <- data_long_foodOnly_factor %>%
  filter(condition == "LO")

##ANOVAs
#HI
hi.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_hi,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(diet),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
hi.kcal.aov$ANOVA
##   Effect DFn DFd        F          p p<.05       ges
## 2   diet   5  50 8.643261 5.8861e-06     * 0.4546631
#LO
lo.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_lo,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(diet),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
lo.kcal.aov$ANOVA
##   Effect DFn DFd        F           p p<.05       ges
## 2   diet   5  50 4.108456 0.003339617     * 0.2793328
#VEH (just so I don't have to scroll back up)
veh.kcal.aov$ANOVA
##   Effect DFn DFd        F            p p<.05      ges
## 2   diet   5  50 8.839876 4.603764e-06     * 0.434701
## grams
# HI dataset
data_grams_long_foodOnly_factor_hi <- data_grams_long_foodOnly_factor %>%
  filter(condition == "HI")

#LO dataset
data_grams_long_foodOnly_factor_lo <- data_grams_long_foodOnly_factor %>%
  filter(condition == "LO")

##ANOVAs
#HI
hi.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_hi,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(diet),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
hi.grams.aov$ANOVA
##   Effect DFn DFd        F            p p<.05       ges
## 2   diet   5  50 10.34611 7.522938e-07     * 0.4985924
#LO
lo.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_lo,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(diet),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
lo.grams.aov$ANOVA
##   Effect DFn DFd        F           p p<.05       ges
## 2   diet   5  50 4.845566 0.001087098     * 0.3126305
#VEH (just so I don't have to scroll back up)
veh.grams.aov$ANOVA
##   Effect DFn DFd        F            p p<.05       ges
## 2   diet   5  50 12.29994 8.511631e-08     * 0.5238583

So we get significant omnibus diet effects for all 3 drug groups (even though the omnibus effect for the low dose is just barely significant). Let’s run pairwise comparisons now:

## HI
#Significance
data_long_foodOnly_factor_hi %>%
  pairwise_t_test(kcal~diet, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
## # A tibble: 15 × 10
##    .y.   group1 group2     n1    n2 statistic    df     p p.adj p.adj.signif
##  * <chr> <chr>  <chr>   <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
##  1 kcal  zero   twenty     11    11   -0.454     10 0.66  1     ns          
##  2 kcal  zero   forty      11    11   -0.362     10 0.725 1     ns          
##  3 kcal  zero   sixty      11    11   -2.10      10 0.062 0.924 ns          
##  4 kcal  zero   eighty     11    11   -4.18      10 0.002 0.028 *           
##  5 kcal  zero   hundred    11    11   -0.0590    10 0.954 1     ns          
##  6 kcal  twenty forty      11    11    0.409     10 0.691 1     ns          
##  7 kcal  twenty sixty      11    11   -1.90      10 0.087 1     ns          
##  8 kcal  twenty eighty     11    11   -3.30      10 0.008 0.119 ns          
##  9 kcal  twenty hundred    11    11    0.452     10 0.661 1     ns          
## 10 kcal  forty  sixty      11    11   -2.31      10 0.044 0.658 ns          
## 11 kcal  forty  eighty     11    11   -3.83      10 0.003 0.05  *           
## 12 kcal  forty  hundred    11    11    0.350     10 0.733 1     ns          
## 13 kcal  sixty  eighty     11    11   -1.73      10 0.115 1     ns          
## 14 kcal  sixty  hundred    11    11    2.14      10 0.058 0.876 ns          
## 15 kcal  eighty hundred    11    11    3.95      10 0.003 0.041 *
#Effect size
data_long_foodOnly_factor_hi %>%
  cohens_d(kcal~diet, paired = TRUE)
## # A tibble: 15 × 7
##    .y.   group1 group2  effsize    n1    n2 magnitude 
##  * <chr> <chr>  <chr>     <dbl> <int> <int> <ord>     
##  1 kcal  zero   twenty  -0.137     11    11 negligible
##  2 kcal  zero   forty   -0.109     11    11 negligible
##  3 kcal  zero   sixty   -0.635     11    11 moderate  
##  4 kcal  zero   eighty  -1.26      11    11 large     
##  5 kcal  zero   hundred -0.0178    11    11 negligible
##  6 kcal  twenty forty    0.123     11    11 negligible
##  7 kcal  twenty sixty   -0.573     11    11 moderate  
##  8 kcal  twenty eighty  -0.996     11    11 large     
##  9 kcal  twenty hundred  0.136     11    11 negligible
## 10 kcal  forty  sixty   -0.695     11    11 moderate  
## 11 kcal  forty  eighty  -1.16      11    11 large     
## 12 kcal  forty  hundred  0.106     11    11 negligible
## 13 kcal  sixty  eighty  -0.520     11    11 moderate  
## 14 kcal  sixty  hundred  0.644     11    11 moderate  
## 15 kcal  eighty hundred  1.19      11    11 large
## LO
#Significance
data_long_foodOnly_factor_lo %>%
  pairwise_t_test(kcal~diet, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
## # A tibble: 15 × 10
##    .y.   group1 group2     n1    n2 statistic    df     p p.adj p.adj.signif
##  * <chr> <chr>  <chr>   <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
##  1 kcal  zero   twenty     11    11    0.737     10 0.478 1     ns          
##  2 kcal  zero   forty      11    11    1.61      10 0.139 1     ns          
##  3 kcal  zero   sixty      11    11   -1.70      10 0.12  1     ns          
##  4 kcal  zero   eighty     11    11   -2.24      10 0.049 0.734 ns          
##  5 kcal  zero   hundred    11    11    0.684     10 0.51  1     ns          
##  6 kcal  twenty forty      11    11    0.187     10 0.855 1     ns          
##  7 kcal  twenty sixty      11    11   -2.21      10 0.051 0.771 ns          
##  8 kcal  twenty eighty     11    11   -3.01      10 0.013 0.198 ns          
##  9 kcal  twenty hundred    11    11    0.0329    10 0.974 1     ns          
## 10 kcal  forty  sixty      11    11   -2.28      10 0.046 0.693 ns          
## 11 kcal  forty  eighty     11    11   -2.90      10 0.016 0.236 ns          
## 12 kcal  forty  hundred    11    11   -0.141     10 0.891 1     ns          
## 13 kcal  sixty  eighty     11    11   -0.0354    10 0.972 1     ns          
## 14 kcal  sixty  hundred    11    11    2.30      10 0.044 0.661 ns          
## 15 kcal  eighty hundred    11    11    2.92      10 0.015 0.229 ns
#Effect size
data_long_foodOnly_factor_lo %>%
  cohens_d(kcal~diet, paired = TRUE)
## # A tibble: 15 × 7
##    .y.   group1 group2   effsize    n1    n2 magnitude 
##  * <chr> <chr>  <chr>      <dbl> <int> <int> <ord>     
##  1 kcal  zero   twenty   0.222      11    11 small     
##  2 kcal  zero   forty    0.484      11    11 small     
##  3 kcal  zero   sixty   -0.513      11    11 moderate  
##  4 kcal  zero   eighty  -0.676      11    11 moderate  
##  5 kcal  zero   hundred  0.206      11    11 small     
##  6 kcal  twenty forty    0.0564     11    11 negligible
##  7 kcal  twenty sixty   -0.667      11    11 moderate  
##  8 kcal  twenty eighty  -0.906      11    11 large     
##  9 kcal  twenty hundred  0.00990    11    11 negligible
## 10 kcal  forty  sixty   -0.686      11    11 moderate  
## 11 kcal  forty  eighty  -0.875      11    11 large     
## 12 kcal  forty  hundred -0.0426     11    11 negligible
## 13 kcal  sixty  eighty  -0.0107     11    11 negligible
## 14 kcal  sixty  hundred  0.694      11    11 moderate  
## 15 kcal  eighty hundred  0.881      11    11 large
## VEH
#Significance
data_long_foodOnly_factor_veh %>%
  pairwise_t_test(kcal~diet, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
## # A tibble: 15 × 10
##    .y.   group1 group2     n1    n2 statistic    df        p p.adj p.adj.signif
##  * <chr> <chr>  <chr>   <int> <int>     <dbl> <dbl>    <dbl> <dbl> <chr>       
##  1 kcal  zero   twenty     11    11    -0.695    10 0.503    1     ns          
##  2 kcal  zero   forty      11    11    -0.704    10 0.498    1     ns          
##  3 kcal  zero   sixty      11    11    -2.10     10 0.062    0.928 ns          
##  4 kcal  zero   eighty     11    11    -4.29     10 0.002    0.024 *           
##  5 kcal  zero   hundred    11    11     2.06     10 0.066    0.996 ns          
##  6 kcal  twenty forty      11    11     0.470    10 0.648    1     ns          
##  7 kcal  twenty sixty      11    11    -1.21     10 0.253    1     ns          
##  8 kcal  twenty eighty     11    11    -3.97     10 0.003    0.04  *           
##  9 kcal  twenty hundred    11    11     1.34     10 0.21     1     ns          
## 10 kcal  forty  sixty      11    11    -1.90     10 0.087    1     ns          
## 11 kcal  forty  eighty     11    11    -4.07     10 0.002    0.034 *           
## 12 kcal  forty  hundred    11    11     2.64     10 0.025    0.37  ns          
## 13 kcal  sixty  eighty     11    11    -1.69     10 0.122    1     ns          
## 14 kcal  sixty  hundred    11    11     2.78     10 0.019    0.29  ns          
## 15 kcal  eighty hundred    11    11     5.26     10 0.000366 0.005 **
#Effect size
data_long_foodOnly_factor_veh %>%
  cohens_d(kcal~diet, paired = TRUE)
## # A tibble: 15 × 7
##    .y.   group1 group2  effsize    n1    n2 magnitude 
##  * <chr> <chr>  <chr>     <dbl> <int> <int> <ord>     
##  1 kcal  zero   twenty   -0.209    11    11 small     
##  2 kcal  zero   forty    -0.212    11    11 small     
##  3 kcal  zero   sixty    -0.634    11    11 moderate  
##  4 kcal  zero   eighty   -1.29     11    11 large     
##  5 kcal  zero   hundred   0.621    11    11 moderate  
##  6 kcal  twenty forty     0.142    11    11 negligible
##  7 kcal  twenty sixty    -0.366    11    11 small     
##  8 kcal  twenty eighty   -1.20     11    11 large     
##  9 kcal  twenty hundred   0.404    11    11 small     
## 10 kcal  forty  sixty    -0.572    11    11 moderate  
## 11 kcal  forty  eighty   -1.23     11    11 large     
## 12 kcal  forty  hundred   0.796    11    11 moderate  
## 13 kcal  sixty  eighty   -0.509    11    11 moderate  
## 14 kcal  sixty  hundred   0.839    11    11 large     
## 15 kcal  eighty hundred   1.59     11    11 large

So without this getting super tedious, it looks like eighty percent takes it across all the two DAMGO conditions as well.

Questions 5 & 6

Keeping in mind that we have a significant main effect of drug on consumption, we can graph consumption of each diet:

data_long_foodOnly %>%
ggplot(mapping = aes(x = factor(condition, level = condition_names), 
                     y = kcal, fill = condition)) + 
         geom_boxplot()

And we can then do pairwise paired t-tests:

## Kcal
# Significance
data_long_foodOnly %>%
  pairwise_t_test(kcal~condition, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
## # A tibble: 6 × 10
##   .y.   group1 group2    n1    n2 statistic    df     p p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
## 1 kcal  HI     LO        66    66     1.40     65 0.167 1     ns          
## 2 kcal  HI     NAL       66    66     2.93     65 0.005 0.028 *           
## 3 kcal  HI     VEH       66    66     2.67     65 0.01  0.058 ns          
## 4 kcal  LO     NAL       66    66     2.04     65 0.046 0.274 ns          
## 5 kcal  LO     VEH       66    66     1.82     65 0.073 0.439 ns          
## 6 kcal  NAL    VEH       66    66    -0.760    65 0.45  1     ns
#Effect size
data_long_foodOnly %>%
  cohens_d(kcal~condition, paired = TRUE)
## # A tibble: 6 × 7
##   .y.   group1 group2 effsize    n1    n2 magnitude 
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>     
## 1 kcal  HI     LO      0.172     66    66 negligible
## 2 kcal  HI     NAL     0.360     66    66 small     
## 3 kcal  HI     VEH     0.328     66    66 small     
## 4 kcal  LO     NAL     0.251     66    66 small     
## 5 kcal  LO     VEH     0.224     66    66 small     
## 6 kcal  NAL    VEH    -0.0935    66    66 negligible
## Grams
data_grams_long_foodOnly %>%
  pairwise_t_test(grams~condition, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
## # A tibble: 6 × 10
##   .y.   group1 group2    n1    n2 statistic    df     p p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
## 1 grams HI     LO        66    66      1.52    65 0.132 0.792 ns          
## 2 grams HI     NAL       66    66      3.03    65 0.004 0.021 *           
## 3 grams HI     VEH       66    66      2.69    65 0.009 0.054 ns          
## 4 grams LO     NAL       66    66      2.19    65 0.032 0.192 ns          
## 5 grams LO     VEH       66    66      1.88    65 0.064 0.386 ns          
## 6 grams NAL    VEH       66    66     -1.00    65 0.319 1     ns
#Effect size
data_grams_long_foodOnly %>%
  cohens_d(grams~condition, paired = TRUE)
## # A tibble: 6 × 7
##   .y.   group1 group2 effsize    n1    n2 magnitude 
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>     
## 1 grams HI     LO       0.188    66    66 negligible
## 2 grams HI     NAL      0.372    66    66 small     
## 3 grams HI     VEH      0.332    66    66 small     
## 4 grams LO     NAL      0.270    66    66 small     
## 5 grams LO     VEH      0.232    66    66 small     
## 6 grams NAL    VEH     -0.124    66    66 negligible

So it looks like the high dose did increase consumption over the vehicle dose, but the low dose didn’t, and naltrexone didn’t significantly lower it. If anything, it looks like drug didn’t lead to huge group level changes, but more just affected the amount of variability across groups.

Additional analysis

We’ve seen so far that drug mildly affects the average amount consumed in general, and that DAMGO exacerbates the preference for the 80% diet while naltrexone flattens it. However, one possibility I want to explore is whether different things are happening for animals that prefer fat vs. those that prefer sugar. This suggests two questions:

7.) Are there differences in baseline preference between animals, and if so, how do we determine those?

8.) Does the pattern of data we obtain differ between fat and sugar preferrers?

In my prospectus, I mention a paper (Gosnell et al., 1990), that describes a procedure for relating baseline preferences to drug effects. I’m going to try to follow their procedure here.

Question 7

As a super simple first pass, I think I can start by filtering the calorie data such that we’re only looking at the vehicle condition and we’re only including the two pure macronutrient conditions:

data_long_foodOnly_pure <- data_long_foodOnly %>%
  filter(condition == "VEH") %>%
  filter(diet == "zero" | diet == "hundred")

data_long_foodOnly_pure %>%
  ggplot(mapping = aes(x = rat,
                       y = kcal,
                       fill = diet)) + 
  geom_col(position = "dodge")

And it’s only at this point that I realize the kcal data probably isn’t the best to be looking at here, because fat is more the twice as caloric as sugar. Let’s repeat the analysis with weight data:

data_grams_long_foodOnly_pure <- data_grams_long_foodOnly %>%
  filter(condition == "VEH") %>%
  filter(diet == "zero" | diet == "hundred")

data_grams_long_foodOnly_pure %>%
  ggplot(mapping = aes(x = rat,
                       y = grams,
                       fill = diet)) + 
  geom_col(position = "dodge")

This still feels like it’s giving us an incomplete picture though. It actually looks like the Gosnell paper had dedicated days from which they detemined baseline preferences. For our data, I can see a few possibilities for which days to use:

1.) Some average of the diet acclimation period. The benefit is that the rats are totally injection naive at this point. Con is that since rats weren’t acclimated to diets in the same order, relative experience in the chambers between fat and sugar diets varies across rats.

2.) An average of the training days where rats were given all 6 diets

3.) The maintenance days in between injections

My gut is telling me that the training days with all 6 diets is the best way to go, so I’m going to start down that road and see what happens.

Let’s load our data:

##Read in data



#kcal
trainingdata_kcal <- read_excel("ID Training Data copy.xlsx", sheet = 32, range = "A20:I31", col_names = column_names[1:9])

trainingdata_kcal_filtered <- trainingdata_kcal
trainingdata_kcal_filtered[trainingdata_kcal_filtered < 0] <- 0

trainingdata_kcal_filtered_long <- trainingdata_kcal_filtered %>%
  pivot_longer(cols = diet_names, names_to = "diet", values_to = "kcal") %>%
  subset(select = -c(h2o, spill))



#grams
trainingdata_grams <- read_excel("ID Training Data copy.xlsx", sheet = 32, range = "A3:I14", col_names = column_names[1:9])

trainingdata_grams_filtered <- trainingdata_grams
trainingdata_grams_filtered[trainingdata_grams_filtered < 0] <- 0

trainingdata_grams_filtered_long <- trainingdata_grams_filtered %>%
  pivot_longer(cols = diet_names, names_to = "diet", values_to = "grams") %>%
  subset(select = -c(h2o, spill))

And filter for pure fat and pure sugar:

trainingdata_kcal_filtered_long_pure <- trainingdata_kcal_filtered_long %>%
  filter(diet == "zero" | diet == "hundred")

trainingdata_grams_filtered_long_pure <- trainingdata_grams_filtered_long %>%
  filter(diet == "zero" | diet == "hundred")

And finally we can visualize

trainingdata_kcal_filtered_long_pure %>%
  ggplot(mapping = aes(x = rat,
                       y = kcal,
                       fill = diet)) + 
  geom_col(position = "dodge")

trainingdata_grams_filtered_long_pure %>%
  ggplot(mapping = aes(x = rat,
                       y = grams,
                       fill = diet)) + 
  geom_col(position = "dodge")

Again, the weight data I think is the best way to go here. Either way, I think a viable strategy here is just to divide rats into fat and sugar preferrers based on the relative amounts they ate here. In that case:

Fat preferrers: ID1, ID12, ID3, ID5, ID9

Sugar preferrers: ID10, ID11, ID2, ID4, ID6, ID7, ID8

Lastly, just to get a sense of each individual rat’s baseline preferences across all diets, I just want to directly plot that as a bar graph:

trainingdata_grams_filtered_long$diet <- factor(trainingdata_grams_filtered_long$diet,
                                       levels = diet_names)

trainingdata_grams_filtered_long$rat <- factor(trainingdata_grams_filtered_long$rat,
                                       levels = rat_names)

trainingdata_grams_filtered_long %>%
  ggplot(mapping = aes(x = diet,
                       y = grams,
                       fill = diet)) + 
  geom_col(position = "dodge") + 
  facet_wrap(~rat)

Based on this, I actually do think a composite measure of baseline preference that accounts of all diets is the best gauge of what the rats actually like. Some of them (like ID12) clearly preferred the highly sugary 80% diet, but didn’t actually eat much of the pure sugar. If we quantified his baseline preferences just from pure sugar, therefore, we’d probably miss out on what his preferences actually were.

Question 8

Part 1

Now let’s try to repeat the analyses from Gosnell et al. In order to repeat their correlation analyses, I need data frames for each non-vehicle diet that contains the following variables:

  • Average baseline fat consumption
  • Average baseline carbohydrate consumption
  • Average baseline for each individual diet
  • A column for each diet reflecting the effect of the drug (i.e., the difference between consumption of the that diet for the current drug dose vs. vehicle)

This is going to be hard, but let’s give it a shot.

We already have a dataset containing each rat’s baseline preferences for each diet, so we can start from there by adding composite fat and sugar variables:

baseline_preferences_grams <- trainingdata_grams %>%
  subset(select = -c(h2o, spill)) %>%
  mutate(
    fat = zero + (twenty*0.8) + (forty*0.6) + (sixty*0.4) + (eighty*0.2),
    sugar = hundred + (eighty*0.8) + (sixty*0.6) + (forty*0.4) + (twenty*0.2)
  )

Now we need to make data frames containing intakes for each drug group:

veh_intakes_grams <- data_grams_filtered %>%
  filter(condition == "VEH") %>%
  subset(select = -c(condition, day)) %>%
  rename(veh_zero = zero,
         veh_twenty = twenty,
         veh_forty = forty,
         veh_sixty = sixty,
         veh_eighty = eighty, 
         veh_hundred = hundred,
         veh_h2o = h2o,
         veh_spill = spill)

hi_intakes_grams <- data_grams_filtered %>%
  filter(condition == "HI") %>%
  subset(select = -c(condition, day)) %>%
  rename(hi_zero = zero,
         hi_twenty = twenty,
         hi_forty = forty,
         hi_sixty = sixty,
         hi_eighty = eighty, 
         hi_hundred = hundred,
         hi_h2o = h2o,
         hi_spill = spill)

lo_intakes_grams <- data_grams_filtered %>%
  filter(condition == "LO") %>%
  subset(select = -c(condition, day)) %>%
  rename(lo_zero = zero,
         lo_twenty = twenty,
         lo_forty = forty,
         lo_sixty = sixty,
         lo_eighty = eighty, 
         lo_hundred = hundred,
         lo_h2o = h2o,
         lo_spill = spill)

nal_intakes_grams <- data_grams_filtered %>%
  filter(condition == "NAL") %>%
  subset(select = -c(condition, day)) %>%
  rename(nal_zero = zero,
         nal_twenty = twenty,
         nal_forty = forty,
         nal_sixty = sixty,
         nal_eighty = eighty, 
         nal_hundred = hundred,
         nal_h2o = h2o,
         nal_spill = spill)

And we need to join the vehicle data frame with that of each other drug dose:

hi_effects_grams <- full_join(veh_intakes_grams,
                              hi_intakes_grams,
                              by = "rat")

lo_effects_grams <- full_join(veh_intakes_grams,
                              lo_intakes_grams,
                              by = "rat")

nal_effects_grams <- full_join(veh_intakes_grams,
                              nal_intakes_grams,
                              by = "rat")

And now we can calculate drug effects by subtracting the intake of each diet on vehicle with the intake of the same diet on drug:

hi_effects_grams <- hi_effects_grams %>%
  mutate(zero_effect = hi_zero - veh_zero,
         twenty_effect = hi_twenty - veh_twenty,
         forty_effect = hi_forty - veh_forty,
         sixty_effect = hi_sixty - veh_sixty,
         eighty_effect = hi_eighty - veh_eighty,
         hundred_effect = hi_hundred - veh_hundred,
         h2o_effect = hi_h2o - veh_h2o,
         spill_effect = hi_spill - veh_spill)

lo_effects_grams <- lo_effects_grams %>%
  mutate(zero_effect = lo_zero - veh_zero,
         twenty_effect = lo_twenty - veh_twenty,
         forty_effect = lo_forty - veh_forty,
         sixty_effect = lo_sixty - veh_sixty,
         eighty_effect = lo_eighty - veh_eighty,
         hundred_effect = lo_hundred - veh_hundred,
         h2o_effect = lo_h2o - veh_h2o,
         spill_effect = lo_spill - veh_spill)

nal_effects_grams <- nal_effects_grams %>%
  mutate(zero_effect = nal_zero - veh_zero,
         twenty_effect = nal_twenty - veh_twenty,
         forty_effect = nal_forty - veh_forty,
         sixty_effect = nal_sixty - veh_sixty,
         eighty_effect = nal_eighty - veh_eighty,
         hundred_effect = nal_hundred - veh_hundred,
         h2o_effect = nal_h2o - veh_h2o,
         spill_effect = nal_spill - veh_spill)

I have a feeling it’ll also make our lives easier if we also join the preference data to each og the effect data frames too:

hi_effects_grams_full <- hi_effects_grams %>%
  left_join(baseline_preferences_grams,
            by = "rat") %>%
  rename(baseline_zero = zero,
         baseline_twenty = twenty,
         baseline_forty = forty,
         baseline_sixty = sixty,
         baseline_eighty = eighty, 
         baseline_hundred = hundred,
         baseline_fat = fat,
         baseline_sugar = sugar)

lo_effects_grams_full <- lo_effects_grams %>%
  left_join(baseline_preferences_grams,
            by = "rat") %>%
  rename(baseline_zero = zero,
         baseline_twenty = twenty,
         baseline_forty = forty,
         baseline_sixty = sixty,
         baseline_eighty = eighty, 
         baseline_hundred = hundred,
         baseline_fat = fat,
         baseline_sugar = sugar)

nal_effects_grams_full <- nal_effects_grams %>%
  left_join(baseline_preferences_grams,
            by = "rat") %>%
  rename(baseline_zero = zero,
         baseline_twenty = twenty,
         baseline_forty = forty,
         baseline_sixty = sixty,
         baseline_eighty = eighty, 
         baseline_hundred = hundred,
         baseline_fat = fat,
         baseline_sugar = sugar)

And now we can actually try to run the correlations! Let’s start with the high condition:

hi.rcorr <- rcorr(as.matrix(hi_effects_grams_full[18:33]))
#hi.rcorr
corrplot(hi.rcorr$r)

Now low:

lo.rcorr <- rcorr(as.matrix(lo_effects_grams_full[18:33]))
#lo.rcorr
corrplot(lo.rcorr$r)

And naltrexone:

nal.rcorr <- rcorr(as.matrix(nal_effects_grams_full[18:33]))
#nal.rcorr
corrplot(nal.rcorr$r)

In the interest of not running a billion significance tests for these correlations, I’m going to draw a subset of the ones I’m most interested in and type them into excel. The file is called “Gosnell results”.

The long and short of it is that this analysis doesn’t suggest that baseline preferences have anything to do with the effects of DAMGO However, I’m going to try to get at this question one more way before I give it up.

Let’s imagine you’re a rat in this study, and you’ve been eating these diets for a while. I think it’s pretty unlikely that the stored representation you have of the diets resembles a thought like “I like fat X amount and sugar Y amount, and will therefore eat from all of the diets in weighted proportion to these preferences.” Instead, I think a much more straightforward way of thinking about it would be something closer to “The sixty percent diet is my favorite, so I’m going to beeline to that one.” So, if DAMGO is really working to exacerbate baseline preferences, what you would expect would be an amplification of that effect. Therefore, there should be a positive correlation between intake of each rat’s favorite diet and the effect of morphine on that diet. Similarly, the correlation should be negative for each rat’s least favorite diet.

Part 2

We’ve already calculated and plotted each rat’s baseline preferences:

trainingdata_grams_filtered_long %>%
  ggplot(mapping = aes(x = diet,
                       y = grams,
                       fill = diet)) + 
  geom_col(position = "dodge") + 
  facet_wrap(~rat,
             ncol = 3)

And it’s at this point that I realize: almost everyone’s favorite diet was the 80%. The DAMGO effect being dependent on baseline preferences is only cool if the rats have different baseline preferences. If their preferences are all the same, then you would expect exactly what we’ve seen – higher doses of DAMGO causing increased consumption of that diet.

This may just be luck of the draw that we got animals who all happened to be 80% preferrers. Alternatively, it may be the case that all rats would prefer 80%. Our dataset doesn’t let us distinguish, but it does reinforce that DAMGO increases consumption of preferred diets.

And just like that, I think we’re done with experiment 2 analysis!

Sike, more plots I want

I want to see each animal’s preferences across testing days too

data_long_foodOnly_factor %>%
  ggplot(mapping = aes(x = diet,
                       y = kcal,
                       fill = diet)) + 
  geom_col(position = "dodge") + 
  facet_wrap(~rat)

data_grams_long_foodOnly_factor %>%
  ggplot(mapping = aes(x = diet,
                       y = grams,
                       fill = diet)) + 
  geom_col(position = "dodge") + 
  facet_wrap(~rat)

I also need to run pairwise comparisons for the one-way ANOVA specifically of the filtered naltrexone data

#Significance
data_grams_long_foodOnly_factor_nal %>%
  pairwise_t_test(grams~diet, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
## # A tibble: 15 × 10
##    .y.   group1 group2     n1    n2 statistic    df     p p.adj p.adj.signif
##  * <chr> <chr>  <chr>   <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
##  1 grams zero   twenty     11    11     0.125    10 0.903 1     ns          
##  2 grams zero   forty      11    11     0.994    10 0.344 1     ns          
##  3 grams zero   sixty      11    11    -1.12     10 0.29  1     ns          
##  4 grams zero   eighty     11    11    -2.64     10 0.025 0.374 ns          
##  5 grams zero   hundred    11    11     1.33     10 0.213 1     ns          
##  6 grams twenty forty      11    11     0.947    10 0.366 1     ns          
##  7 grams twenty sixty      11    11    -1.19     10 0.26  1     ns          
##  8 grams twenty eighty     11    11    -2.89     10 0.016 0.243 ns          
##  9 grams twenty hundred    11    11     1.07     10 0.31  1     ns          
## 10 grams forty  sixty      11    11    -1.64     10 0.132 1     ns          
## 11 grams forty  eighty     11    11    -4.11     10 0.002 0.032 *           
## 12 grams forty  hundred    11    11     1.03     10 0.329 1     ns          
## 13 grams sixty  eighty     11    11    -0.666    10 0.521 1     ns          
## 14 grams sixty  hundred    11    11     1.76     10 0.109 1     ns          
## 15 grams eighty hundred    11    11     4.42     10 0.001 0.019 *

I need to reformat some data for easier use in SPSS

data_kcal_wide <- data_long_foodOnly %>%
  subset(select = -day) %>%
  pivot_wider(names_from = c(condition, diet), values_from = kcal)
#write_xlsx(data_kcal_wide, "/Users/lawilson1999/Library/Mobile Documents/com~apple~CloudDocs/Psychology MA/Second Year (2022-2023)/Thesis/Thesis analyses (ID)\\data_kcal_wide.xlsx")

data_grams_wide <- data_grams_long_foodOnly %>%
  subset(select = -day) %>%
  pivot_wider(names_from = c(condition, diet), values_from = grams)
#write_xlsx(data_kcal_wide, "/Users/lawilson1999/Library/Mobile Documents/com~apple~CloudDocs/Psychology MA/Second Year (2022-2023)/Thesis/Thesis analyses (ID)\\data_grams_wide.xlsx")

I want to run additional tests comparing the effect of drug dose across diet groups for kcals:

###Filter out diet groups
#kcals
# 0% dataset
data_long_foodOnly_factor_zero <- data_long_foodOnly_factor %>%
  filter(diet == "zero")

#20% dataset
data_long_foodOnly_factor_twenty <- data_long_foodOnly_factor %>%
  filter(diet == "twenty")

#40%
data_long_foodOnly_factor_forty <- data_long_foodOnly_factor %>%
  filter(diet == "forty")

#60%
data_long_foodOnly_factor_sixty <- data_long_foodOnly_factor %>%
  filter(diet == "sixty")

#80%
data_long_foodOnly_factor_eighty <- data_long_foodOnly_factor %>%
  filter(diet == "eighty")

#100%
data_long_foodOnly_factor_hundred <- data_long_foodOnly_factor %>%
  filter(diet == "hundred")

##ANOVAs
#0%
zero.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_zero,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
zero.kcal.aov$ANOVA
##      Effect DFn DFd        F         p p<.05       ges
## 2 condition   3  30 1.436595 0.2516532       0.0223938
#20%
twenty.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_twenty,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
twenty.kcal.aov$ANOVA
##      Effect DFn DFd         F        p p<.05        ges
## 2 condition   3  30 0.9382934 0.434407       0.01115226
#40%
forty.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_forty,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
forty.kcal.aov$ANOVA
##      Effect DFn DFd        F        p p<.05        ges
## 2 condition   3  30 1.086801 0.369725       0.06477213
#60%
sixty.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_sixty,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
sixty.kcal.aov$ANOVA
##      Effect DFn DFd        F          p p<.05        ges
## 2 condition   3  30 2.580396 0.07197252       0.08632641
#80%
eighty.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_eighty,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
eighty.kcal.aov$ANOVA
##      Effect DFn DFd        F            p p<.05       ges
## 2 condition   3  30 9.205299 0.0001795536     * 0.3107644
pairwiseEighty <- data_long_foodOnly_factor_eighty %>%
  pairwise_t_test(kcal~condition, paired = TRUE,
                  p.adjust.method = "bonferroni",
                  estimate = TRUE)
pairwiseEighty
## # A tibble: 6 × 10
##   .y.   group1 group2    n1    n2 statistic    df     p p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <dbl> <chr>       
## 1 kcal  HI     LO        11    11     2.08     10 0.065 0.388 ns          
## 2 kcal  HI     VEH       11    11     2.28     10 0.046 0.277 ns          
## 3 kcal  HI     NAL       11    11     3.14     10 0.01  0.062 ns          
## 4 kcal  LO     VEH       11    11     0.554    10 0.592 1     ns          
## 5 kcal  LO     NAL       11    11     1.52     10 0.159 0.954 ns          
## 6 kcal  VEH    NAL       11    11     1.99     10 0.075 0.448 ns
#100%
hundred.kcal.aov <- ezANOVA(data = data_long_foodOnly_factor_hundred,
                    dv = .(kcal),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
hundred.kcal.aov$ANOVA
##      Effect DFn DFd        F         p p<.05       ges
## 2 condition   3  30 1.683584 0.1915543       0.1022568

And for grams

###Filter out diet groups
#kcals
# 0% dataset
data_grams_long_foodOnly_factor_zero <- data_grams_long_foodOnly_factor %>%
  filter(diet == "zero")

#20% dataset
data_grams_long_foodOnly_factor_twenty <- data_grams_long_foodOnly_factor %>%
  filter(diet == "twenty")

#40%
data_grams_long_foodOnly_factor_forty <- data_grams_long_foodOnly_factor %>%
  filter(diet == "forty")

#60%
data_grams_long_foodOnly_factor_sixty <- data_grams_long_foodOnly_factor %>%
  filter(diet == "sixty")

#80%
data_grams_long_foodOnly_factor_eighty <- data_grams_long_foodOnly_factor %>%
  filter(diet == "eighty")

#100%
data_grams_long_foodOnly_factor_hundred <- data_grams_long_foodOnly_factor %>%
  filter(diet == "hundred")

##ANOVAs
#0%
zero.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_zero,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
zero.grams.aov$ANOVA
##      Effect DFn DFd        F         p p<.05       ges
## 2 condition   3  30 1.436595 0.2516532       0.0223938
#20%
twenty.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_twenty,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
twenty.grams.aov$ANOVA
##      Effect DFn DFd         F        p p<.05        ges
## 2 condition   3  30 0.9382934 0.434407       0.01115226
#40%
forty.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_forty,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
forty.grams.aov$ANOVA
##      Effect DFn DFd        F        p p<.05        ges
## 2 condition   3  30 1.086801 0.369725       0.06477213
#60%
sixty.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_sixty,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
sixty.grams.aov$ANOVA
##      Effect DFn DFd        F          p p<.05        ges
## 2 condition   3  30 2.580396 0.07197252       0.08632641
#80%
eighty.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_eighty,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
eighty.grams.aov$ANOVA
##      Effect DFn DFd        F            p p<.05       ges
## 2 condition   3  30 9.205299 0.0001795536     * 0.3107644
#100%
hundred.grams.aov <- ezANOVA(data = data_grams_long_foodOnly_factor_hundred,
                    dv = .(grams),
                    wid = .(rat),
                    within = .(condition),
                    return_aov = TRUE,
                    type = 2)
## Warning: You have removed one or more Ss from the analysis. Refactoring "rat"
## for ANOVA.
hundred.grams.aov$ANOVA
##      Effect DFn DFd        F         p p<.05       ges
## 2 condition   3  30 1.683584 0.1915543       0.1022568

Scatters showing correlation between preference and intake:

hi_effects_grams_full %>%
  ggplot(mapping = aes(x = baseline_eighty,
                       y = eighty_effect)) +
  geom_point() +
  geom_smooth(method = "lm") + 
  labs(title = "High Dose (0.25µg/0.5µL/side)",
       x = "Baseline 80% Sugar Intake (grams)",
       y = "Effect of DAMGO on 80% Sugar Intake (grams)") +
  theme_classic(base_size = 16)
## `geom_smooth()` using formula = 'y ~ x'

lo_effects_grams_full %>%
  ggplot(mapping = aes(x = baseline_eighty,
                       y = eighty_effect)) +
  geom_point() +
  geom_smooth(method = "lm", se = 0.95) + 
  labs(title = "Low Dose (0.025µg/0.5µL/side)",
       x = "Baseline 80% Sugar Intake (grams)",
       y = "Effect of DAMGO on 80% Sugar Intake (grams)") +
  theme_classic(base_size = 16)
## `geom_smooth()` using formula = 'y ~ x'

Let’s try to get these scatterplots on the same plot:

hi_effects_eighty <- hi_effects_grams_full %>%
  subset(select = c(rat, baseline_eighty, eighty_effect)) %>%
  mutate(dose = "HI")

  

lo_effects_eighty <- lo_effects_grams_full %>%
  subset(select = c(rat, baseline_eighty, eighty_effect)) %>%
  mutate(dose = "LO")

damgo_effects <- rbind(hi_effects_eighty, lo_effects_eighty)

damgo_effects %>%
  ggplot(mapping = aes(x = baseline_eighty,
                       y = eighty_effect,
                       color = dose)) +
  geom_point() +
  geom_smooth(method = "lm") + 
  labs(
       x = "Baseline 80% Sugar Intake (grams)",
       y = "Effect of DAMGO on 80% Sugar Intake (grams)") +
  theme_classic(base_size = 13) +
  scale_color_discrete(name = "DAMGO Dose", labels = c("High (0.25 µg/0.5µL/side)", "Low (0.025µg/0.5µL/side"))
## `geom_smooth()` using formula = 'y ~ x'